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We develop a consistent semiclassical method to calculate the probability of topo¬ 
logical soliton-antisoliton pair production in collisions of elementary particles. In our 
method one adds an auxiliary external field pulling the soliton and antisoliton in the op¬ 
posite directions. This transforms the original scattering process into a Schwinger pair 
creation of the solitons induced by the particle collision. One describes the Schwinger 
process semiclassically and recovers the original scattering probability in the limit of 
vanishing external field. We illustrate the method in (1 -|- 1)-dimensional scalar field 
model where the suppression exponents of soliton-antisoliton production in the multi¬ 
particle and two-particle collisions are computed numerically. 
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1 Introduction 

An intriguing possibility of observing nonperturbative phenomena in particle collisions and in 
the early Universe triggered development of powerful semiclassical methods for description of 
false vacuum decay [1, 2, 3], instanton-like transitions [4, 5] and multiparticle production [6, 
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Figure 1: (a), (b) Production of the soliton-antisoliton (SA) pair in the two-particle collision, 
(c) A tree-level diagram with 2Ns particles in the hnal state. 

7]. Although these processes occur with exponentially small rates [8, 9, 10, 11, 12, 13, 14] in 
weakly interacting theories, they may be important for baryogenesis [15, 5], stability of the 
Higgs vacuum [16], or violation of the baryon and lepton numbers at future colliders [17]. 

We consider nonperturbative transitions of a new type: classically forbidden creation of 
topological soliton-antisoliton pairs in W-particle collisions, see the sketch of the N = 2 
process in Figs. la,b. At weak coupling the probabilities of these transitions are expected 
to be exponentially suppressed. Indeed, the solitons can be crudely regarded [18] as bound 
states of Ns oc l/g"^ particles, where g 1 is the coupling constant. The cross section 
of their pair production in the collision of few particles involves the multiparticle factor [6] 
Vfew ~ 5'^^® (2Ws)! ~ exp(—const/^f^) due to roughly {2Ns)\ tree-level diagrams exemplihed 
in Fig. Ic. This argument, however hand-waving^, suggests a general expression for the 
probability 

Vn{E) ~ (1.1) 

of the inclusive process N particles —)■ solitons + particles; here Fn [E) is the multiparticle 
suppression exponent. The suppression should disappear at sufficiently large number of col¬ 
liding particles N > 2Ns when creation of the solitons proceeds classically. Presently the 
form (1.1) is supported by unitarity arguments [9] and recent calculation [19] in a quantum 
mechanical model describing the soliton moduli space. However, no reliable held theoretical 
method for evaluating the exponent F]s{E) is known. Below we propose a general semiclas- 
sical method of this kind which is applicable, in particular, to the processes with N = 2 
initial particles. 

^In particular, loop corrections are shown [6] to be large at g^Ns ^ 1- 
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Figure 2: Soliton-antisoliton pair in the external field £. 

Recently several dynamical mechanisms for enhancing the rate of soliton-antisoliton pro¬ 
duction in particle collisions were proposed [20]. They are supported by classical arguments 
which cannot be directly extended to quantum level. Our method for computing the expo¬ 
nent in Eq. (1.1) will be valuable for tests of these mechanisms and their phenomenological 
applications. 

So far the processes in Figs. la,b eluded semiclassical treatment. The reason can be 
traced back to the attraction between the topological soliton and antisoliton [ 21 ]: taken 
at rest, they accelerate towards each other and annihilate^ into many particles. Thus, the 
potential barrier between the soliton-antisoliton pair and multiparticle states is absent. As 
a consequence, classically forbidden transitions from particles to solitons cannot be directly 
described by the powerful semiclassical methods developed for tunneling processes. 

In this paper which extends Ref. [22] we solve the above difficulty by relating production of 
the solitons from particles to their Schwinger pair creation in external field. We illustrate our 
method with explicit numerical calculations and present details of the numerical techniques. 

The above relation is established in the following way. By definition the topological 
soliton and antisoliton can be equipped with the topological charges ±g; we denote by 
the respective topological current, = 0. We introduce an external U{1) field 

coupled to this current^, 

AR = - y , (1.2) 

where d is the number of spacetime dimensions. We consider the field = (—0) with 
the constant “electric” component S which pulls the soliton and antisoliton in the opposite 
directions, see Fig. 2. Then the solitons are Schwinger pair produced like the opposite charges 
in the ordinary electric field. In our method one starts from the well-known semiclassical 
solution describing spontaneous creation of the soliton pair at 7 ^ 0. Then one adds N 
particles with energy E to the initial state and obtains solutions for the collision-induced 

^We do not consider integrable models where annihilation of the soliton and antisoliton may be forbidden 
by conservation laws. 

^The interaction in Eq. (1.2) does not need to be consistent at the quantum level. It is introduced as an 
auxiliary semiclassical tool and will be switched off in the end of the calculations. 
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Schwinger processes, cf. [23]. Finally, one raises the particle energy above the kinematic 
threshold 2Ms of soliton-antisoliton production, where Ms is the soliton mass, and takes 
the limit of vanishing external held —>■ 0. In this way one arrives at the semiclassical 
solutions describing creation of the solitons from particles. The semiclassical suppression 
exponent Fjq{E) is computed using these solutions. 

To be concrete, imagine that our method is applied to pair production of ’t Hooft- 
Polyakov monopoles [24] in particle collisions. The modihcation term (1.2) in this case cou¬ 
ples gauge-invariant magnetic current [25] to a small external held £. One can show [26] 
that the latter is equivalent to external magnetic held. Then the monopole-antimonopole 
pairs are created in the Schwinger process in accordance with the electric-magnetic duality. 
One therefore obtains the semiclassical solutions relevant for the monopole-antimonopole 
production in particle collisions from the solutions [26, 27] describing their creation in the 
magnetic held. 

The main technical difficulty of our method is related to the semiclassical description 
of the collision-induced Schwinger processes. Below we use the approach of [28, 11, 12] 
which involves a family of complex semiclassical solutions satisfying a certain boundary value 
problem. The suppression exponent Fn{E) is calculated as a functional on these solutions. 
In the end of the calculation we send T —)■ 0. 

An obstacle to the semiclassical method appears in the phenomenologically interesting 
case of = 2 colliding particles: the two-particle initial state cannot be described semiclassi- 
cally. We overcome this obstacle by appealing to the Rubakov-Son-Tinyakov conjecture® [28] 
which states that the multiparticle exponent F^{E) does not depend on the initial particle 
number N at semiclassically small values of the latter i.e. a.t N ^ l/g"^, see [30, 31, 32] for 
con&mations. This means that the two-particle exponent F 2 {E) coincides with Fsf{E) in 
the limit g'^N —)■ 0. Calculating semiclassically' the latter exponent and taking the limit, we 
obtain F 2 {E). 

We illustrate the above semiclassical method in the (1-1-1)-dimensional scalar held model. 


S = 


1 


dt dx 


I (^4'f 


VW 


(1.3) 


where the coupling constant® g is small, = {t,x), and the scalar potential R(0) has two 
degenerate vacua 0 = see Fig. 3a, solid line. The model (1.3) possesses a pair of kink-like 

®An alternative strategy involving singular solutions is suggested in [29, 10, 7]. 

^Recall that the semiclassical description is valid at 1, which can hold even in the limiting region 
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Figure 3: (a) Scalar potential V{(j)) (solid line) and its modification (dashed line), (b) Soliton 
(S) and antisoliton (A) solutions, (c) Potential energy of their interaction. 


static solutions interpolating between the vacua: topological soliton (S) and antisoliton (A) 
shown in Fig. 3b. Below we consider their pair production in the A^-particle collisions. 

Following the general strategy outlined above, we introduce an external field 
= {—Sx, 0) coupled to the topological current® where £ is negative in 

accordance with Fig. 2 and is the antisymmetric symbol in two dimensions. Integrating 
by parts, we rewrite the modification term (1.2) as 

AS =—£ j (j)dtdx. (1.4) 

This interaction can be absorbed in the redefinition of the scalar potential V ^ V + g‘^£4>- 
After that 0_ and 0+ become false and true vacua, V{(j)-) > 14(0+) . One concludes that the 
Schwinger pair production of the kink-like solitons in the external field £ is equivalent [2] to 
the spontaneous false vacuum decay via the nucleation of true vacuum bubbles. The latter 
is a celebrated tunneling process with well-known semiclassical description [3]. 

In what follows we modify the scalar potential 14(0) giving negative energy density 
14(0+) = —6p to the true vacuum, 6p oc \£\. We calculate the rate of false vacuum de¬ 
cay accompanied by the A^-particle collision [33, 10, 11, 14] and remove the modification 

—)■ 0 in the end of the calculation. 

Let us visualize the potential barrier between the soliton-antisoliton pair and multiparti¬ 
cle sector of the modified system. Kink-like soliton and antisoliton at distance I attract each 
other with Yukawa force Fatt oc e“'"+0 where m\ = 14"(0+), cf. [21]. Besides, they are pushed 
apart by the pressure 5p. This leads to the potential energy of the static soliton-antisoliton 

g'^N < 1 . 

®Field rescaling (f> ^ g(f> brings this constant in front of the interaction terms in the action. 

®The related topological charge is / Jodx = 0(x = -|-oo) — 0(x = —oo). 
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Figure 4: (a) Regimes of transitions in the [E, N) plane of initial data, (b) Numerical results 
for the suppression exponent F]^{E) of soliton-antisoliton production in the model (1.3). 


pair 

[/(/) ~—I7oe“”*+* —/dp + const (1.5) 

shown in Fig. 3c. The barrier U{1) separates the multiparticle sector of the theory from the 
static soliton-antisoliton pairs with / > l^r, where Icr ~ |log(const ■ dp)|/m+ is the critical 
distance corresponding to the barrier top. The pairs with / < Icr annihilate and therefore 
should be attributed to the multiparticle sector. Unstable static solution (j)cb{x) describing the 
soliton and antisoliton at a distance I = Icr is the famous critical bubble [1, 2, 3]; its energy 
Ecb determines the height of the potential barrier between the solitons and multiparticle 
states. As 6p —)■ 0, the critical bubble turns into infinitely separated soliton and antisoliton 
at rest: Icr +oo and Ecb 2 M 5 . In this limit the potential barrier is hidden below the 
kinematic threshold E = 2Ms of soliton pair production. 

False vacuum decay at U = = 0 (no initial particles) is described by the bounce 

solution [3]. In the main body of the paper we numerically^^ hnd similar semiclassical 
solutions 0s(/, x) for the decay accompanied by the collision of N particles at energy E. We 
observe several regimes of transitions summarized in Fig. 4a. The shaded region E < Nm_, 
where m_ is the particle mass in the false vacuum, is kinematically forbidden. In the opposite 
case of high energies and N > 4/ (region III in Fig. 4a) transitions proceed classically and 
we obtain F/v(U) = 0. The latter classical regime is investigated in [34], see also [35]. 

At smaller multiplicities the false vacuum decay is classically forbidden, Fn{E) > 0. We 
find that the properties of the respective semiclassical solutions (j)s(t, x) are qualitatively 

use the scalar potential shown in Fig. 3a; expression for V{(j)) will be given in the main text. 
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different at energies below and above the height Ecb of the potential barrier (regions I and II in 
Fig. 4a). Solutions with E < E^b resemble the bounce: they describe formation of large true 
vacuum bubbles and imply inhnite suppression^^ Ep^{E) oc 5p~^ in the limit 5p —?■ 0. High- 
energy solutions from the region II in Fig. 4a have smaller spatial extent. They approach 
the critical bubble (j)cbix) plus outgoing waves as t —)■ +oo. Thus, at inhnitesimally small 
6p they describe creation of widely separated soliton-antisoliton pairs at rest. We explicitly 
demonstrate that the respective value of the suppression exponent En{E) approaches a hnite 
value as 6p —)■ 0. 

Our numerical results for the suppression exponent of producing the soliton pair from 
particles aX 6p = 0 are shown in Fig. 4b. The two-particle exponent E 2 {E) is obtained by 
evaluating the additional limit g‘^N —)■ 0. One observes that at high energies E]sf{E) is hnite 
and (almost) energy-independent which is the expected behavior [9, 10, 14]. 

The paper is organized as follows. We introduce the semiclassical method in Sec. 2 and 
numerical technique in Sec. 3. Numerical solutions and results for the suppression exponent 
are described in Sec. 4. We conclude in Sec. 5. Technical details are presented in Appendices. 


2 The semiclassical boundary value problem 


Consider false vacuum decay induced by a collision of N particles at energy E. The inclusive 
probability of this process equals 


V, 


N 




{f\U{tf, U)\z- E,N) 


,-FM{E)/g^ 


( 2 . 1 ) 


where U is the quantum evolution operator and the inhnite time limits tij —>■ =Fcx) are 
assumed. The sum in Eq. (2.1) runs over all V-particle initial states with energy E in the 
false vacuum and hnal states |/) containing a bubble of true vacuum. In the approximate 
equality we introduced the suppression exponent En{E), c.f. Eq (1.1). 

At small the action (1.3) is parametrically large and any path integral with exp(iS') in 
the integrand can be evaluated in the saddle-point approximation. In Appendix A we review 
the semiclassical method for calculating the probability (2.1): introduce a path integral for 
Vn{E) and derive equations for its (generically complex) saddle-point conhguration (j)s(t, x), 
see Refs. [28, 37, 5] for details. The leading semiclassical exponent exp(—F/v/^f^) is given by 
the value of the integrand at 0 = 0^. 

^^This property can be deduced from the thin-wall approximation [33, 36] which predicts large suppression 
at low energies but breaks down a,t E > 2Ms. 





Figure 5: (a) Contour in complex time for the semiclassical boundary value problem. Crossed 
circles represent singularities of the solution, (b) Schematic form of the semiclassical solution 
(j)sit,x). Pink/gray region corresponds to 0 ~ 0+. 

Let us describe the saddle-point equations for x) which will be solved in the next 
Sections. This complex conhguration should extremize the action (1.3) i.e. satisfy the clas¬ 
sical held equation 

- diet), + = 0 , ( 2 . 2 a) 

where V = dV/d(j). We will consider Eq. (2.2a) along the complex-time contour in Fig. 5a 
corresponding to evolution of particles prior to collision (part AB of the contour), tunneling 
via formation of the true vacuum bubble (part BC) and expansion of the bubble to inhnite 
size (part CD). The expected structure of the semiclassical solution along the contour is 
visualized in Fig. 5b. The duration T of Euclidean evolution is a parameter of the solution. 

The held equation is supplemented with the boundary conditions related to the initial 
and hnal states of the process. Namely, in the asymptotic future the solution should be real, 

Im 0 s, —)■ 0 as t ^+oo , ( 2 . 2 b) 

due to the inclusive hnal state in Eq. (2.1). In the asymptotic past i.e. a.t t = iT +1' and 
t' —)■ — oo, the conhguration (p, should evolve linearly in the false vacuum, 

OO 

0s —>• + [ —= COS (kx) (fk 0 —>■ —oo , (2.2c) 

J Jnujk ^ ' 

0 
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where ul = k'^ + and we expect that the relevant semiclassical solutions are x —)■ —x 
symmetric like the soliton pair in Fig. 3b. The waves in Eq. (2.2c) represent free on- 
shell particles in the initial state of the process. The initial saddle-point condition relates 
negative- and positive-frequency components of these waves, 

fk = e~^gk , (2.2d) 


where 6 is the second parameter of the solution. In the special case 6 = 0 Eq. (2.2d) implies 
real-valued (pg at t' —?• —oo. The initial state in this case is inclusive i.e. optimized with 
respect to the multiplicity N at hxed energy E, cf. Eq. (2.2b). The semiclassical solutions 
at 6^ = 0 are called periodic instantons [38], we will consider them in the Sec. 4. At 0 > 0 
the solutions are complex at t' —)■ —oo. Moreover, at 6* —)■ +oo Eq. (2.2d) reduces to the 
Feynman boundary condition /^ = 0. The initial state in this limit is indistinguishable from 
the vacuum and therefore contains semiclassically small number of particles, N <C 

Equations (2.2) form the T/9 boundary value problem [28] which will be used to hnd the 
semiclassical solutions (ps- The parameters T and 9 of the solutions are Lagrange multipliers 
conjugate to the energy E and initial multiplicity respectively. The latter quantities are 
given by the standard expressions 

OO OO 

g'^E = 2 J dkujkfkgl, g^N = 2 j dk fkgl- (2.3) 

0 0 

In what follows we hnd (ps for all possible values of T and 9 and compute {E, N) by Eqs. (2.3). 

The suppression exponent in Eq. (2.1) is evaluated as a functional on the semiclassical 
solution <ps(t, x), 


fjv(£) = 2/1111 S|^J-y(2ET + Af«) + 2Im dx(^,-4,.)d4, , 


(2.4) 


where the last two terms represent contributions from the nontrivial initial state of the 
process. Note that the semiclassical parameter g does not enter the boundary value prob¬ 
lem (2.2). Thus, the exponent (2.4) depends on the combinations g^E and g‘^N rather than 
individually on E and N, cf. Eqs. (2.3). This feature is in agreement with the Rubakov- 
Son-Tinyakov conjecture [28] discussed in the Introduction: the semiclassical exponent does 
not depend on N ai N <^1/g'^ the limit 


E 2 {E) = lim EjsiiE) 

g'^N^O 


(2.5) 
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exists. Then F 2 {E) is the suppression exponent of the two-particle processes. 

We remark that the method of Lagrange multipliers implies the following Legendre trans¬ 
forms for T and 6, 


2/T = -dEF^iE) , g^e = -dNF^iE) , (2.6) 

see Appendix A for derivation. Below we will keep in mind that T and 6 are proportional 
to the partial derivatives of Fn{E). 

Let us comment on the delusively simple hnal state of the process. First, recall that 
the conhgurations 0s(t, x) should describe false vacuum decay i.e. contain the bubble of 
true vacuum — expanding or critical — at t —)■ -|-oo (pink/gray region in Fig. 5b). This 
property is nontrivial and will be used as a selection rule for the relevant semiclassical 
solutions. Second, one might think that the semiclassical solutions are real at the real 
time axis: starting from real 0^ and dt4>s in the asymptotic future and solving the classical 
held equation backwards in time, one arrives at real 0s (t, x) at t E M. This simple logic is, 
however, not applicable if the solution approaches the critical bubble (pcbix) in the asymptotic 
future [39]. Indeed, the latter is unstable and contains exponentially growing and decaying 
modes (50±(f,x) cc in the spectrum of linear perturbations. At large positive times 

one obtains 0s ~ 0c6 + A_(50_ -|- real waves, and the coefficient A_ is complex in general. 
Then the overall solution, although complex-valued at the real time axis, satishes Eq. (2.2b) 
asymptotically. We will demonstrate that the semiclassical solutions with E > Ect approach 
the critical bubble in the inhnite future. In this case the asymptotic reality (2.2b) cannot 
be imposed at hnite real t; we overcome this difficulty using the methods of [39, 40, 32]. 

To summarize, we arrived at the practical method for evaluating the suppression exponent 
of soliton pair production in particle collisions. One starts by describing the induced false 
vacuum decay at dp > 0: hnds a family of the solutions 0s(t, x) to Eqs. (2.2), relates their 
parameters (T, 9) to {E, N) by Eqs. (2.3) and computes the suppression exponent i0v(E) 
using Eq. (2.4). The exponent of the soliton-antisoliton production is recovered in the limit 
dp ^ 0 above the kinematic threshold E > 2Ms- Besides, in the limit g‘^N —)• 0 one obtains 
the exponent F 2 {E) of the two-particle processes, see Eq. (2.5). 
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3 Numerical methods 


3.1 Choosing the potential 


Nonlinear semiclassical equations (2.2) do not admit analytic treatment and we solve them 
numerically. To this end we specify the scalar potential in dimensionless units, 


v{<p) = ^{<p+ir 


l-vW 


- 1 


where W{u)=e (m + + m®) (3.1) 


and a = 0.4. This function has a double-well form with false and true vacua at 0 = 0_ = — 1 
and (;/) = (/)+> 0, respectively. We have V{cj)_) = 0 and £x the energy density V{(f)+) = —6p 
of the true vacuum by tuning the parameter v. In particular, at n ~ 0.75 the vacua are 
degenerate; the respective scalar masses are m_ ~ 1, m_|_ 7.6. Below we start from the 

solutions at (5p > 0 (larger v), then send 6p —)■ 0. Function (3.1) is shown in Fig. 3a at (5p = 0 
and 6p = 1 (solid and dashed lines, respectively). The soliton and antisoliton conhgurations 
in the case of degenerate vacua are demonstrated in Fig. 3b. 

Let us explain the choice (3.1) of the scalar potential. First, we do not consider the 
standard 0^ theory because evolution of the kink-antikink pairs there is known to exhibit 
chaotic behavior [41]. That chaos is related to the fact that the spectrum of linear perturba¬ 
tions around the 0^ kink contains two localized modes representing its spatial translations 
and periodic vibrations of its form. The modes accumulate energy during the kink-antikink 
evolution which is therefore described by two^^ collective coordinates. This evolution is 
chaotic like the majority of two-dimensional mechanical motions. Irregular dynamics is dif¬ 
ficult [42, 43, 44, 45] for the semiclassical methods, we avoid it by changing the form of 
the scalar potential. We checked that the model (3.1) is not chaotic^^, i.e. there is a single 
localized mode in the spectrum of linear perturbations around the soliton. 

Second, the potential (3.1) is almost quadratic at 0 < 0 and essentially nonlinear at 
positive 0. This property is convenient for numerical methods. Indeed, the initial con¬ 
ditions (2.2c), (2. 2d) should be imposed in the asymptotic past where 0s(t, x) < 0 de¬ 
scribes wave packets linearly evolving in the false vacuum (diagonal lines in the left part of 
Fig. 5b). In the model (3.1) the wave packets remain free almost up to their collision point, 
and the initial condition can be imposed at relatively small negative times (realistically, at 
Keti ~ —(6 7)/m_). Long nonlinear evolution in other models can be costly for numerical 


^^The other two coordinates representing center-of-mass motion and P-odd vibrations of the kink- 
antikink pair decouple due to momentum conservation and x —>■ —x symmetry. 

^^Another type of irregular behavior [46] appearing at < m_ is also not met in our model. 
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Figure 6 : Rectangular lattice superimposed with the schematic form of the semiclassical 
solution j. The time sites tj cover the contour ABCD in Fig. 5a. Solution at x < 0 is 
reconstructed using the x —>■ —x symmetry. 

methods. Note that this problem with slow linearization is specihc to (l+l)-dimensional 
systems and should be absent in higher dimensions. 


3.2 Discretization 


We discretize the semiclassical equations (2.2) using the rectangular (W + 3) x lattice 
in Fig. 6 . The sites tj and Xi of the lattice cover the contour ABCD in Fig. 5 and the 
interval 0 < x < respectively. The solution at x < 0 can be restored from the reflection 
symmetry x) = —x). Our lattice is as close to uniform as possible. Namely, 

Xi = i ■ Ax with integer 0 < i < A^^, — 1 and spacing Ax = L^j{N^ — 1). The time spacings 
Atj = — tj, where —1 < j < Nt + 1, are constant at the real and imaginary parts of 

the contour. Typically, we choose \Atj\ Ax. This property will be useful for formulating 
the boundary conditions. Our unknowns (j)j^i = (j)s(tj, Xj) are the (complex) values of the 
semiclassical solution at the lattice sites. 

We replace the derivatives in the action (1.3) with the second-order hnite-difference 
expressions 


dt(j)\ 


i+i/2p 


. -)■ 


At. 


dx(j)\ 


iP+i/2 


-)■ 


bj, i+l 4^j, 
Ax 


(3.2) 


taking values at the links; below we mark all link quantities with half-integer indices. We 
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also trade the time and space integrals for the second-order trapezoidal snms, e.g. 

. Nt iVt + 1 

fj+i/2 - E Atj fj , (3.3) 

i=-i i=-i 

where the hrst and second expressions are nsed for the link and site qnantities, respectively. 
The last snm involves the “averaged” spacings Aij = {Atj + Atj_i)/2 and Ai_i = At_i/2, 

At^t-^i = AtNt/‘2- 

Snbstitntions (3.2), (3.3) tnrn the action into a fnnction S'(0) of {Nt + 3) x complex 
variables For completeness we present this fnnction in Appendix B. The lattice held 
eqnations are then obtained by extremizing S'(0) with respect to the held valnes at the 
interior sites of the contonr ABCD: 


dS/d(t)j^i = 0 , 


0 < j < Nt and all i . 


(3.4) 


Note that the eqnations at the last spatial site i = N^ — 1 correspond to the Nenmann 
bonndary condition dx(j) = 0 imposed at x = L^. Note also that Eqs. (3.4) can be derived 
by writing the lattice path integral for the probability and evalnating it in the saddle-point 
approximation, like in the continnons approach of Appendix A. 

The hnal bonndary condition (2.2b) implies that the held is real at the last two time 
sites, 

= Im07Vt+i,i = 0 for all i , (3.5) 

cf. Eq. (3.2). In Sec. 2 we remarked that this condition cannot be nsed for solutions ap¬ 
proaching the critical bnbble at f —>■ -1-cxd. That case will be considered separately. 

Discretization of the initial conditions (2.2c), (2. 2d) is far less trivial. Let ns hrst simplify 
the discnssion and consider the system in discrete space {xt} and continnons time t. Recall 
that the semiclassical evolntion is linear in the beginning of the process t ~ t_i. We therefore 
introdnce a deviation of the held from the false vacnnm 'ipi{t) = [(j){t, Xi) — Axt and 
leave only qnadratic terms in the action. We arrive at the system of N^ conpled oscillators. 


(7^5(2) = 


. N^-1 'j 

I i=0 i,i'=0 J 


(3.6) 


where the real symmetric matrix ifj j/ is a discretized version of the operator {—d^ + m^), 
see Appendix B for the explicit form. Linear evolntion in Eq. (3.6) is solved by decomposing 
'?/’i(t) in real-valned eigenvectors of H with eigenvalnes o;^. 


Mt) = 


N, 




(n) 


^ ^/2uJr 


n 




(3.7) 
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Here we introduced the integration constants /„, g-n and shifted the time by iT to keep 
contact with Eq. (2.2c). By direct inspection one observes that the the index n, vectors 
and constants /„, g* are the lattice analogs of k, cos{kx) and fk, gl in Eq. (2.2c), 
respectively. Thus, the discrete version of the initial condition (2.2d) is fn = Q~^gn- The 
quantities entering this condition are extracted from as 


fn 





1 






(3.8) 


where we used Eq. (3.7) assuming normalization ^ = ^n,n' and omitting irrel¬ 

evant phase factor. To obtain the lattice form of the initial condition, we discretize the time 
derivatives in Eqs. (3.8). The simplest way is to perform the substitutions (3.2), (3.3) in the 
quadratic action (3.6), and then dehne dt'ipi at the very Erst time site as 


dti^i 


t-1 


g^ 

2 d%l^i{t_i) 


Af_i 


+ Af_i */?/>*/(f_i) , 

i' 


(3.9) 


cf. the derivation of the initial conditions in Appendix A. 

We see that fn and gn in Eqs. (3.8), (3.9) are the linear combinations of V’i(^-i), 'fiifo) 
and their complex conjugates. This means that the initial condition /„ = gn can be 
explicitly rewritten as a set of 2Nx real equations on T = 'ififto)}, 


Mr ■ Re T M/ ■ Im T = 0 . 


(3.10) 


Constant 2Nx x 2Nx matrices Mr, Mi can be deduced from Eqs. (3.8), see Appendix B. In the 
numerical code we explicitly compute^^ the eigenvectors and eigenvalues of Hi ii, construct 
matrices Mr, Mi, and add Eqs. (3.10) to the system of discrete equations. 

Discretization leaves us with {Nt -|- 1) x complex held equations (3.4) and real 
boundary conditions (3.5), (3.10). This is precisely the required number of equations to 
determine {Nt -|- 3) x complex unknowns 0^,^. The semiclassical equations will be solved 
in the subsequent Sections. 

Finally, let us introduce lattice expressions for the energy, initial particle number and 
suppression exponent. The full nonlinear energy 

g^E= [ dx[{dtcl)f + {d^cj)f + 2V{(l))] (3.11) 

Jo 

^"^This is done once for each lattice. 
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is directly discretized using the substitutions (3.2), (3.3). We monitor its conservation to 
control the discretization errors. The initial energy and particle number (2.3) in the discrete 
case have the form, 

g^E = 2j2^^fn9:, g^N = 2j2fng:, (3.12) 

n=l n=l 

see Appendix B. The suppression exponent Fn{E) is computed using the discrete action 
S{(j)j^i) and quantum numbers (3.12) in Eq. (2.4); the last integral in this expression is 
discretized in the standard way. 

3.3 Fixing the time translations 

The continuous semiclassical equations (2.2) preserve time^^ shifts. Namely, (j)s{t — to, x) 
with real to satishes the equations if (f>sit, x) does. This feature, if left unnoticed, leads to 
a numerical artifact. Indeed, since the time translations are explicitly broken in the lattice 
system, the position to of the numerical solution is out of direct control: it is hxed by the 
discretization and hnite-volume effects. Depending on the lattice parameters, it may become 
arbitrary large and cause divergence of the numerical procedure. 

To cure the above artifact, we notice [12] that due to the continuous time-shift sym¬ 
metry one of the lattice equations starts to depend on the others at Af —>■ 0, t_i —)■ —oo. 
Indeed, in this limit the total energy of the solution is conserved. It is real due to the final 
boundary conditions (3.5) and, on the other hand, proportional to the sum Yhn ^nfnQn 
Eq. (3.12). Now, consider the initial conditions /„ = e~^gn which imply, in particular, 
relations arg/„ = axggn- One of the the latter is redundant because the sum Unfngn is 
automatically real. We can exclude the redundant equation without affecting the solution. 

Following this line of reasoning, we drop the phase of the initial condition /„q = e~^gno 
at a given n = no relating only the absolute values 

\f^^\ = e-%r.o\ • (3.13) 

To keep the number of lattice equations equal to the number of unknowns, we add one 
condition hxing the time translation invariance of the solution. Provided the set of the 
lattice equations is solved, the arguments of /„p and gn^ will be automatically equal up 
to discretization errors. In realistic numerical calculations we use Eq. (3.13) for a highly 

^^The spatial translations are fixed by cc —> —x symmetry imposed on the semiclassical solutions. 


16 



populated mode with^® no = 2 or 3; the equality of phases is then satished with accuracy 
better than |arg/„g — arg^f^^l < 10“^ for all our solutions. 

Let us now discuss the additional equation hxing the time translation invariance. It is 
convenient to use different conditions at energies below and above the barrier height Ecb- At 
low energies we impose the relation [11] 

dxw{x)Redt(j)\^^ = 0 with w{x) = (3.14) 

at the point C of the time contour in Fig. 5a. Equation (3.14) places the turning point 
dt4> = 0 of the solution, if there is one, at a given time t = tc- Even if the solution does 
not have turning points, Eq. (3.14) keeps the singularities of the solution (crossed circles in 
Fig. 5a) away from the time contour. 

At high energies the structure of the semiclassical solutions changes, and we use different 
constraint for the time translations. Namely, we hx the center-of-mass coordinate R of the 
wave packet at the start of the process, see Fig. 6, 

Lx 

dx (r — i?) Rep£;(a;) ^ = 0 , (3.15) 

where pe is the energy density entering Eq. (3.11): g‘^E = 2 dx pe{x). In realistic 
calculations we take R ~ O.TLa, to guarantee that the initial wave packet is inside the lattice 
range at t = f_i. 

Constraints (3.14), (3.15) are discretized in the standard way using the substitutions 
(3.2) and (3.3). We repeat that exchanging one complex equation fn^ = Q~^gno for the two 
real — Eq. (3.13) and one of Eqs. (3.14), (3.15) — one keeps the number of equations equal 
to the number of unknowns. 

3.4 Solving the equations 

Let us select the appropriate values of the lattice parameters. To this end we regard the 
schematic form of the semiclassical solution 0^ j in Fig. 6. 

First of all, the site numbers Nt, should be large enough to ensure small discretization 
errors O(Af^), O(Ax^). Note, however, that the characteristic frequencies of the solutions are 
estimated as Uk ^ E/N, see Eqs. (3.12). They become higher at larger energies and smaller 
multiplicities. In these regions the solutions are sharp and require hne lattice resolution. On 

^®The respective physical momentum is « 7r(no — RIL^. 
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the other hand, large Nt and are costly for numerical algorithms. In practical calculations 
we use lattices NtX = 3200 x 150 and 11000 x 500 below and above Ecb, respectively. We 
keep track of the discretization errors and stop^^ obtaining solutions whenever the relative 
accuracies become worse than 1%. 

Next, we want to keep the nonlinear part of the solution at x. Re t ~ 0. Then the initial 
time Re ti = should be large negative to ensure linear evolution at the start of the process. 
We find that^® Ref_i = -(6-^8) is large enough: in this case the relative contributions of 
nonlinear interactions at t ~ t_i are smaller than 10“^. Besides, the spatial volume |x| < 
should encompass the entire solution. Since the initial wave packets propagate inside the 
lightcone (diagonal lines in Fig. 6), we require > |Ret_i|. In practice it is convenient to 
fix^® = 7, then choose t_i to keep the wave packets at |x| < L^. We checked that our 
results are insensitive to the value of at the relative level 10“^. Finally, we explained in 
Sec. 2 that the solutions with E < E^b are real at the real time axis. In this case we reduce 
the part CD of the time contour to two sites where the condition (3.5) is imposed. 

At E > Ecb the solution becomes real-valued only asymptotically at t —)■ +oo, and we are 
obliged to extend the contour to large positive times. We choose ~ |Ret_i|. Then the 
waves emitted in the interaction region remain within the lattice range, see Fig. 6. 

We performed several tests of the numerical procedure, see Appendix C for details. The 
overall conclusion is that the linearization and hnite-volume effects cause relative errors 
smaller than 10“^, while the relative discretization accuracies remain below 1%. 

We proceed by describing the numerical procedure to solve the set of algebraic lattice 
equations. Our choice of the numerical method is limited because the field equation (2.2a) is 
of hyperbolic and elliptic types at the Minkowski and Euclidean parts of the time contour, 
respectively. In addition, the semiclassical solution 0s(t, x) is complex and satisfies the 
boundary rather than initial conditions. The most convenient numerical technique in this 
case [11, 12] is based on multidimensional Newton-Raphson method. To simplify notations, 
let us denote the complete set of lattice field equations (3.4) and boundary conditions (3.5), 
(3.10), (3.13), (3.15) by Ei^k{.4>) = 0. In the Newton-Raphson method one starts from 
the approximation 0^°] to the solution and repeatedly refines it by finding the correction 

^^See [14] for the study of the high-energy region. 

^®Recall that in our units to_ « 1. 

At low energies we use = 8 15 due to larger sizes of the respective true vacuum bubbles. 
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6(f) = (ps — 0^°^- The latter is obtained from the set of linear eqnations 




^ d(pj,i 

or ’ 


4>m 


^4^j, i 


0 


(3.16) 


in the backgronnd of After obtaining the correction, one redefines the approximation 
—)■ + 6(f> and repeats the procednre. The iterations stop once Sep becomes smaller 

than the predetermined nnmerical error. 

The Newton-Raphson method converges qnadratically [47]. Typically, the acceptable 
precision 6(p < 10“^° is achieved in 5 — 6 interactions. However, this method is extremely 
sensitive to the very hrst approximation (p^^'^\ even slightly incorrect choice of the latter may 
canse divergence of the iterative procednre. 

We therefore nse a carefnl strategy for finding the nnmerical solntions. First, we obtain 
the “simpler” solntion at some particnlar valnes of the parameters (T, 9) = (Tq, 9q). The 
obvions candidate is the bonnee [3]. It satisfies the semiclassical bonndary valne problem 
(2.2) at Tq = + 00 , 00 = 0 and, on the other hand, can be compnted using half-analytic 
methods. In the next Section we will discuss an extended one-parametric family of “simpler” 
solutions. Second, we apply the Newton-Raphson method to hnd the solution at (T, 6) = 
(Tq + AT, 6q + A6) using the one at (Tq, 0o) as the first approximation. At sufficiently small 
AT, A9 the method has to converge. Third, starting from the new solution, we change 
(T, 6) by a small step, again, and obtain yet another solution. We repeat this procedure 
until a complete two-parametric family of the numerical solutions is found. 

We finish this Section with a remark that the Newton-Raphson method requires nu¬ 
merical solution of the sparse linear system (3.16). This is the most time-consuming part 
of the numerical procedure. We perform it using two alternative algorithms described in 
Appendix D. Our “fast” algorithm arrives at the solution in CPU time tcpu oc W • oper¬ 
ations but has poor stability properties. We find that it works for the high-energy solutions 
but accumulates round-off errors at T < Ecb, see explanation in Appendix D. In the latter 
region we use stable (yet slower) algorithm involving tcpu oc W • operations. Both our 
algorithms are highly parallelizable and implemented at the multiprocessor cluster. 
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4 Numerical results 


4.1 Periodic instantons 

We start by considering the periodic instantons [38] — semiclassical solutions at 0 = 0 and 
arbitrary T. The boundary conditions (2.2b), (2.2c) and (2.2d) in this case imply reality 
of (ps{t, x) in the beginning and at the end of the process i.e. along the Minkowski parts 
AB and CD of the time contour in Fig. 5a. Further insight is gained if we assume that the 
solutions have turning points at the corners B and C of the contour, 

9t0s(O, x) = dt(f)s{iT, x) = 0 . (4.1) 

Then the semiclassical evolution along the Euclidean part BC also proceeds with real-valued 
0s(t, x). The solutions satisfying Eq. (4.1) are 2T-periodic in Euclidean time. 

Physical arguments [38] suggest that all relevant periodic instanton solutions satisfy the 
Ansatz (4.1). Indeed, consider a somewhat different process: transition between the sectors 
of the false and true vacua at temperature Its rate F(/d) is obtained by integrating the 
multiparticle probability (2.1) with the Boltzmann exponent, 

F(/?) ~ ; (4.2) 

where the prefactors are ignored. The integral (4.2) is saturated near the saddle point 

P + dEF^/g^ = /? - 2T = 0 , = -g^O = 0 , (4.3) 

where we used the Legendre transforms (2.6) in the hrst equalities. One sees that the saddle- 
point parameters 0 = 0 and T = /d/2 in Eq. (4.3) correspond to a periodic instanton^^. Using 
Eq. (2.4) at real 0^, we obtain the thermal rate 

Y r-u g-2Im5[(^q 

suppressed by the Euclidean action of this solution. One hnds it natural that real Euclidean 
solutions with period f3 = 2T describe thermal transitions. 

In our study the very hrst periodic instanton is obtained from the critical bubble 4>cb{x) — 
the static true vacuum bubble at the verge of collapse. In Fig. 7 we show ((>cb{x) at 6p = 0.4. 
Since the critical bubble is unstable, the spectrum of linear perturbations in its background 

such solution exists. At high temperatures /? < 27r/|w_|, where oj- is defined below, the integral (4.2) 
is saturated near the barrier top E = Ecb- 
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Figure 7: Critical bubble (j)cb{x) and the spectrum of linear perturbations in its back¬ 
ground; 6p = 0.4 and L^; = 7. 

contains an exponentially growing mode 6(j)-{x) with eigenfrequency cai < 0. This means 
that the conhguration 

x) = (j)cb{x) + A_ 5(f)_{x) cosh(|a;_|t) (4.5) 

solves the held equation at small A_. The approximate solution (4.5) is periodic in Euclidean 
time with turning points at t = 0 and t = 7rz/|a;_|. It satishes the boundary condition (4.5) 
and therefore represents the periodic instanton with T = 7r/|a;_|. 

To compute numerically the approximate conhguration (4.5), we solve the static held 
equation and hnd the critical bubble (j)cb{x)- The spectrum of its linear perturbations is 
given by the matrix-diagonalization routine of Sec. 3.2, see the inset in Fig. 7. Picking up 
the negative mode with < 0 and choosing A_ = 0.3, we construct the approximate 
solution (4.5). 

Now, we are ready to hnd the entire family of the periodic instantons. We numerically 
solve the held equation along the Euclidean part BC of the time contour with the boundary 
conditions (4.1). We use the Newton-Raphson method described in the previous Section. 
The very hrst solution (Fig. 8c) is obtained at T = 7r/|a;_| -|- AT, where AT is small. In 
this case the conhguration (4.5) serves as a zeroth-order approximation for the numerical 
procedure. Next, we increase the parameter T in small steps hnding one solution at a time. 
The starting conhguration at each step is provided by the last known solution. Examples of 
the periodic instantons are plotted in Figs. 8a-c. Their energies E{T) are given by Eq. (3.11) 
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Figure 8: (a)—(c) Three-dimensional plots of the periodic instantons x) G M at different 
T; these solutions have 6^ = 0. Only Euclidean parts t G [iT, 0] of the solutions are shown. 
Red/dark and blue/light regions correspond to > 0 and 0^ < 0, respectively. Lower 
panel: Parameters T{E) of the periodic instantons. Circles with letters correspond to the 
solutions (a)-(c). 


(squares in the lower panel of Fig. 8). 

One sees a clear distinction between the low- and high-energy solutions in Figs. 8a and 8c. 
While the latter resembles the critical bubble, the former is nearly rotationally-invariant and 
has large Euclidean period T. At E —)■ 0 the solutions approach the bounce thus describing 
spontaneous decay of the false vacuum. The periodic instantons are absent in the opposite 
region E > E^t- 
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Figure 9: Euclidean action of the periodic instantons at different 6p (dashed lines) and their 
limit —)■ 0 (empty points) versus the thin-wall result (4.6) (solid line). 


Now, consider the limit —)■ 0 of the solutions with 6* = 0. In Appendix E we remind 
that the sizes of the true vacuum bubbles become inhnite at small 5p justifying the celebrated 
thin-wall approximation [3, 33, 10, 36, 48]. The respective solutions 0s(t, x) can be obtained 
analytically. Their action equals 


2Im S'[0s] 


6p 


arcsin 


TSp 


, T6p I fT6py 
WMs) 


at T < g‘^Ms/5p (4.6) 


and stays constant at larger periods, see Appendix E and [36]. Note that the typical values of 
T and 2Im S in Eq. (4.6) are proportional tol/6p implying inhnite suppression F]^{E) —)■ +cx) 
at Sp —?• 0. This comes with no surprise, since transitions between the degenerate vacua are 
energetically forbidden below the threshold E = 2Ms of soliton pair creation. On the other 
hand, at hnite T and —)■ 0 the rate (4.4), (4.6) reduces to the Boltzmann probability 


T ~ 


g-2M,/3 




of Ending the soliton pair in the thermal ensemble at temperature I3~^ = (2T)“^, cf. [49]. 

In Fig. 9 we compare the actions of the periodic instantons at different 6p (dashed lines) 
with Eq. (4.6) (solid line). One observes nice agreement which becomes almost perfect if one 
extrapolates the numerical results to = 0 (empty points in Fig. 9). 
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Figure 10: Solutions at E < E^b- (a) The periodic instanton with (T, 9) = (5.2, 0), 
g‘^{E, N) ^ (5.0, 3.8). (b) Real part of the solution Re0s(t, x) at (T, 6) = (5.0, 0.76) and 
{g'^E, g'^N) k. (5.0, 2.0). Only the central regions of moderately small |t| and |x| < 5 are 
shown. The points A, B, C, D of the time contour are marked above each solution. Color 
represents arg^^. 

4.2 Solutions below Ecb 

Solving the held equation backwards and forwards in time, we continue the periodic instan- 
tons to the parts AB and CD of the time contour. We thus obtain the complete solutions, see 
the one in Fig. 10a. At large negative times they describe wave packets in the false vacuum 
which collide at f ~ iT producing expanding bubbles of the true vacuum. We compute the 
in-state quantum numbers {E, N) of the solutions by Eqs. (2.3). The line of the periodic 
instantons is shown with hlled squares in the (i7, N) plane of Fig. 11. 

Starting from the known solutions at 6 = 0, we hnd the ones with positive 6. To this end 
we increase the value of 6 in small steps keeping T = const. At each step we numerically 
solve^^ the boundary value problem (2.2). An example of the solution 'Re(j)s{t, x) with 0 > 0 
is shown in Fig 10b. It is complex and contains sharp wave packets at large negative time, 
see the color representing arg^^. 

Evaluating the quantum numbers {E, N) of the solutions by Eqs. (2.3), we mark them 
with points in the initial data plane of Fig. 11. One sees that the numerical data cover the 
region E < Ecb and N > lA/g"^-. at small N the solutions become sharp and require better 

^^Recall that the lattice analogs (3.5), (3.14) of the final boundary conditions are imposed at the point C 
of the time contour. After solving the equations, we continue each solution to the part CD of the contour. 
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Figure 11: Parameters {E, N) of the numerical solutions at < Ecb- Empty and filled 
points are situated along the lines T = const and 9 = const, respectively; the value of the 
constant parameter is written near each line in black for T and in green/gray for 6. Circles 
with letters represent parameters of the solutions in Fig. 10. Filled squares are the periodic 
instantons (6* = 0). 


lattice resolution. On the other hand, the high-energy region will be explored in the next 
Section. The suppression exponent En{E) is computed using Eq. (2.4). 

The thin-wall arguments of the previous Section suggest that the exponent is inversely 
proportional to 6p at small values of the latter. Thus, its Laurent series expansion starts 
with the singular term, 

E^{E) = + F^^oiE) + 0{5p) . (4.7) 

dp 

In Appendix E we deduce 


EM,-i{E)=2g^Ml 


arccos 




(4.8) 


considering dynamics of the thin-wall bubbles, see also [36]. The function (4.8) is shown by 
solid line in Fig. 12. It does not depend on N and decreases with energy reaching zero at 
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Figure 12: The thin-wall prediction (4.8) (solid line) and numerical results for the two- 
particle exponents F 2 {E) at E < 2Ms and different dp (dashed lines). Points are obtained 
by linearly extrapolating the numerical data to hp = 0. 

E = 2Ms- At higher energies Eq. (4.8) does not give any sensible result: in Appendix E we 
argue that the thin-wall approximation breaks down in that region. Presumably, this means 
that the limit —)■ 0 of the suppression exponent exists at E > 2Ms- 

Let us compare the numerical results for En{E) with the thin-wall expression (4.8). 
Since the periodic instantons with 9 = 0 have been already studied, we consider the opposite 
case 6 —?• +oo or g‘^N —)• 0. In this limit Fn{E) coincides with the exponent F 2 {E) of 
transitions from the two-particle initial states, see the conjecture (2.5). We extrapolate the 
numerical data for Fn {E) to g‘^N = 0 and obtain the dashed lines in Fig. 12; the details 
of this extrapolation will be discussed in Sec. 4.6. One observes that the numerical graphs 
approach Eq. (4.8) at smaller 6p and coincide with it after the additional extrapolation to 
6p = 0 (empty points). 

4.3 Going to £” > Ecb 

Since the thin-wall bubbles do not describe transitions at E > Ecb, one expects to find new 
properties of the semiclassical solutions in that region. In realistic calculations we observe 
somewhat different effect: at energies above^^ Ecb our numerical method either diverges 
or produces unphysical “reflected” solutions exemplified in Fig. 13. The latter satisfy the 
semiclassical equations but approach the false vacuum at f —)■ +cxd. They apparently cannot 

^^More precisely, this happens at E > Ei{N) « Ecb- We ignore small difference between these two 
thresholds in what follows. 
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Figure 13: The unphysical “reflected” solution at > Ecb] (T, 6) = (1.866, 0.3) and 
g‘^{E, N) ^ (6.11, 4.1). 

describe transitions between the vacua. One immediately identifles [39] the root of the 
problems: a condition forcing the solutions to interpolate between the two vacua is not 
present in the semiclassical boundary value problem (2.2). As a consequence, the numerical 
procedure can produce unphysical solutions even if the ones with correct properties exist. 

The problem of fixing the qualitative behavior of the saddle-point configurations is rather 
general. We solve it using the e-regularization technique of Refs. [39, 40, 32]. To this end we 
recall that in our numerical method the solutions with shrinking bubbles are continuously 
obtained from the physical ones by decreasing the parameter T below some value T^{9). The 
“boundary” solutions at T = E{9) are very special: their true vacuum bubbles neither shrink 
nor expand but survive to f —>■ +cxd. The main idea of the e-regularization is to exclude 
all “boundary” configurations from the set of accessible semiclassical solutions. Once this 
is achieved, one cannot obtain solutions with shrinking bubbles from the physical ones by 
continuous deformations. Then aX E > E^b we will And solutions with correct properties (if 
they exist). 

With this logic in mind, we add [39] a small imaginary term to the classical action, 

S' —>■ S'e = S'[(/)] + ieTintlfj)] , (4-9) 

where the modification parameter e and functional Tint are positive-definite. Importantly, 
we require special properties of Tj„t[0]: it should take finite values on any configuration 
interpolating between the vacua and diverge if x) contains a static finite-size bubble 
in the Anal state. Then the latter “boundary” configurations with static bubbles represent 
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singularities of the modified action S'e —>■ +ioo and cannot coincide with its extrema. To 
the contrary, the semiclassical solutions in the modihed system extremize and do not 
belong to the class of “boundary” conhgurations. We conclude that the modihed solutions 
cannot change qualitative properties in the course of continuous deformations; hnding their 
continuous family at £' > Ect and sending e —)■ +0, one arrives at correct high-energy 
solutions. 

Let us construct the functional Tj„t[0] with the desired properties. We choose 

Tint[0] = f{x)Wint ’ (4-10) 

where f{x) = e~^ and Wint{u) = a = 0.4. The function f{x) restricts the spatial 

integral to the central region |x| < a/, where cx/ = 0.4 in our numerical calculations. At the 
same time, Wint is vanishingly small at 0 ~ and takes positive values between the vacua, 
0_ < 0 < 0_|_. Accordingly, the time integral in Eq. (4.10) diverges if the conhguration 
x) contains a static hnite-size bubble at t —)■ +cxd. However, if the conhguration is 
physical and the bubble expands, the value of the held at |x| < <7/ tends to 0+ at large times 
and the f-integral converges. We conclude that the functional (4.10) discriminates between 
the “boundary” and physical conhgurations. Note that the modihcation (4.9), (4.10) simply 
deforms the scalar potential. 


H ^ K(</>, x) = H(0) 


ief{x)Wint 


(^) 


( 4 . 11 ) 


introducing minimal changes of the numerical code. 

We remark that the regularization (4.9) is pretty general: it was successfully applied in 
quantum mechanics [39, 40, 32, 45, 50], held theory [12] and gravity [51]. Besides, it can be 
justihed by adding a constraint to the path integral with the Faddeev-Popov trick [40, 32], 
cf. [52]. 

Following the above recipe, we pick up some physical solution at F < Ecb, introduce 
regularization (4.11) with e = 5 ■ 10“^ and hnd the modihed solution with the same T and 
9. After that we decrease T obtaining the set of modihed solutions at F > Ecb, see Fig. 14a 
and cf. Fig. 13. As expected, all these solutions contain expanding bubbles at f —)■ +cxd 
and therefore describe transitions between the vacua. We hnally decrease the modihcation 
parameter to inhnitesimally small values e ;< 10“^, Fig. 14b, and continue to explore the 
plane of initial data by changing the parameters T and 6 in small steps. Eventually, we 
obtain all possible solutions at F > Ecb, see Fig. 15. 





Figure 14: Modified solutions at 6 = 0.3, E 6.11/g‘^ > E^b and different values of 
the regularization parameter: (a) e = 5 ■ 10“^, (b) e = 10“^. The other parameters are 
(a) T = 0.75, g^N ^ 4.14, (b) T = 0.7, g'^N ^ 4.13. 

Figure 16 displays two solutions at E > Ecb] their quantum numbers {E, N) are marked 
by circles with letters in Fig. 15. The initial wave packets in these solutions are composed 
of high-frequency modes, their true vacuum bubbles are small, cf. Fig. 10. Besides, the 
respective periods T of Euclidean evolutions are short. This allows us to use the faster 
numerical algorithm of Appendix D. Recall that the solutions become singular in the limits 
of small N and high E] for our best lattice this bounds the accessible region to > l/^f^ 
and E < U/g'^. 

We stress that the transition to the region E > E^b is smooth at e > 0. In particular, the 
parameters N and T of the modihed solutions change smoothly along the lines 9 = const, see 
Fig. 17. However, this crossover becomes sharper at smaller e suggesting continuous rather 
than differentiable changes in the e = 0 solutions across Ecb- 

Now, consider the limit e —)■ +0 of the semiclassical solutions. Figure 14 displays two 
solutions with E > Ecb at different values of the regularization parameter e. The true 
vacuum bubble in the solution with smaller e remains static for a notably longer period 
before it starts to expand. This suggests that the original solution with e = 0 contains a 
static hnite-size bubble in the hnal state {t —)■ +cxd). The latter property becomes apparent 
once we plot the t = const sections of the solution Re0s(t, x) with inhnitesimally small e, 
see Fig. 18a, dashed lines. At large times this solution approaches the critical bubble (solid 
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Figure 15: A complete set of numerical solutions in the {E, N) plane; 6p = 0.4, notations 
of Fig. 11 are used. For simplicity the results a.t E > Ecb are obtained at nonzero e < 10“^ 
using the smaller lattice Nt x = 3200 x 150. Solutions with the highest energies g'^E > 11 
and smallest multiplicities g‘^N <1.6 are not shown: they need better lattice resolution. 
Circles with letters are the solutions in Figs. 16a, b. 

line): 

(psit, x) —)■ (j)cb{.x) + outgoing waves as f —>■ +oo . (4-12) 

Importantly, we hnd that the asymptotics (4.12) holds for all original (e = 0) solutions above 
Ecb, not just the ones at the boundary of this region. Thus, all high-energy solutions are 
unstable and cannot be obtained by direct numerical method. 

The unusual behavior of the solutions at £' > Ecb suggests that the respective processes 
of false vacuum decay proceed in two stages. First, the critical bubble is created with 
exponentially small probability; the energy excess E — Ecb is dropped in the form of the 
outgoing waves. Second, the critical bubble, being classically unstable, decays producing a 
growing bubble of true vacuum with probability of order one. Similar tunneling mechanisms 
appear in multidimensional quantum mechanics [44, 39, 40, 32, 50] and other models of held 
theory [12, 13]. 

Finally, consider the limit e —?• +0 of the quantum numbers (i?, N) and exponent E^ 
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Figure 16: Solutions at > E^b, e = 5 ■ 10“^ and 6p = 0.4. Their parameters are indicated 
by circles with letters in Fig. 15. Namely, the solution (a) has (T, 6) = (1.8 ■ 10“^, 0.1) and 
g'^{E, N) ^ (8.0, 3.8), in the case (b) (T, 0) = (8 ■ 10“^, 0.6) and g‘^{E, N) ^ (10.0, 2.1). 




Figure 17: Transition to the high-energy region E > E^b at e > 0. The graphs show the 
parameters N (left panel) and T (right panel) of the modified solutions as functions of energy 
E along the line 9 = 0.3. Large circles mark the solutions in Figs. 14a,b. 


at hxed T and 9. We find that these quantities are regular functions of e: the data points 
in Fig. 18b are well fitted with the quadratic polynomials (dashed lines). We quadratically 
extrapolate them to e = 0 and obtain the exponent En{E) of false vacuum decay at high 
energies. The numerical errors related to this extrapolation procedure are negligibly small. 
The suppression exponent is plotted at g'^N = 2 and Sp = 0.4 in Fig. 19. 
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Figure 18: (a) Fixed-time sections of the numerical solution Re0s(t, x) with E > Ecb and 

inhnitesimally small e. Numbers near the graphs are the values of Ref. Solid line shows 
the critical bubble ((>cb{x). The parameters of the solution are (T, 9) = (10“^, 0.1) and 
g‘^{E, N) ^ (8.6, 3.8); dp = 0.4. (b) The limit e —)■ +0 of the suppression exponent F/v and 
energy E at (T, 6) = (0.01, 0.6), g‘^{E, N) ^ (9.5, 2.1) and dp = 0.4. Dashed lines are the 
quadratic hts of the data points. 



Figure 19: The suppression exponent Eisf{E) of collision-induced false vacuum decay at 
high energies and g‘^N = 2] dp = 0.4. The graph is obtained in the limit e —)■ +0 using the 
largest lattice Nt x = 11000 x 500. 
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Figure 20: Parameters {E, N) corresponding to the classically allowed transitions at 6p = 
0.4 (dots) and their lower boundary Nmin{E) (dashed line). The solid line is the same 
boundary extracted from the semiclassical results. 

4.4 Classical over—barrier transitions 

An important test of our semiclassical method involves classically allowed decay of the false 
vacuum hlled with many particles. These processes proceed with unsuppressed probabilities, 
Fn {E) = 0, their initial quantum numbers {E > Ecb, N) occupy the upper right corner in 
Fig. 15. Importantly, the minimal particle number N = Nmin{E) required for such transitions 
can be obtained in two ways: by studying the real classical solutions and by hnding the region 
E]sf{E) = 0 from the classically forbidden side. Comparing the two results, one checks the 
semiclassical method at high energies. 

We studied the classically allowed decay of the multiparticle states in the false vacuum 
in [34] using the stochastic sampling technique of [35]. Namely, we constructed many sets 
of random Cauchy data {0(ti, x), dt(f>(ti, a:)} in the false vacuum and obtained a classical 
solution for each set. We selected the solutions arriving to the true vacuum, i.e. those con¬ 
taining expanding bubbles at f —>■ -l-oo. Finally, we calculated the initial quantum numbers 
{E, N), Eqs. (2.3), for the selected solutions and indicated them with green/gray dots in 
Figs. 15 and 20. We arrived at the region in the {E, N) plane corresponding to the classi¬ 
cally allowed decay of the false vacuum; its lower boundary (dashed line in Fig 20) gives the 
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Figure 21: (a) The semiclassical solution at 5p = 0.03; its parameters (i?, N) are the 
same as in Fig. 16b. (b) Linear extrapolation of the suppression exponent to Sp = 0 at 
g^(E, N) = (9, 2). 


minimal initial multiplicity Nmin{E) of classical over-barrier transitions. 

On the other hand, the function Nmin{E) can be obtained using the semiclassical results 
of the previous Section. One notes that any real solution to the classical held equations 
satishes the boundary value problem (2.2) at T = 6 = 9 and gives En{E) = 0. Thus, 
at T, 6* —)■ +0 the semiclassical solutions become real and their quantum numbers {E, N) 
approach the boundary N = Nmin{E) of classical over-barrier transitions. Besides, since 
En{E) vanishes at the latter boundary. 


dN„ 


dE 


dN 

~dE 


Fn(E)=0 


lim 

T, 0-5-+O 


2T 

T 


( 4 . 13 ) 


where the Legendre transforms (2.6) were used in the last equality. Relation (4.13) implies 
that the limit T, 0 —)■ 0 should be performed at '6 = 6/T = const, where -6 parametrizes 
the curve Nmin{E). An exemplary semiclassical solution with small T and 9 is plotted in 
Fig. 16a, its quantum numbers are marked by the circle in Fig. 15. 

In Fig. 20 we plot the boundary Nmin{,E) extracted from the semiclassical results (solid 
line). It almost coincides with that obtained from the real classical solutions. This supports 
our semiclassical method in the high-energy region E > E^b- 
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Figure 22: Suppression exponent F^{E) at g'^N = 2 and different 6p. Solid line is obtained 


by linearly extrapolating the data to Sp = 0. 

4.5 Soliton—antisoliton production: > 0 

So far we have considered the collision-induced decay of the false vacuum. In this Section 
we send 6p —)■ +0 and extract the exponent Fn{E) of the soliton-antisoliton pair production 
in particle collisions. 

Keeping E > E^b and e > 0, we decrease the energy density of the true vacuum to 
6p 0.02 0.06. We arrive at the semiclassical solutions exemplihed in Fig. 21a. Their 
true vacuum bubbles expand at small constant velocities without any apparent acceleration. 
Sending, in addition, e —)■ +0, we obtain the solutions arriving at the static critical bubble 
(pcbix). Recall that the spacial size of (pcb is logarithmically large at small 6p. 

As expected, we hnd that the semiclassical exponent Fn (E) is not very sensitive to in 
the high-energy region E > Ecb- Moreover, at small 6p its dependence is linear, see Fig. 21b. 
This means that the singular term in the “thin-wall” expansion (4.7) is absent and the limit 
—)■ 0 of the exponent exists. Linearly extrapolating Fj^f{E) to 6p = 0 at hxed E and 
N (dashed line in Fig. 21b), we obtain the exponent of soliton-antisoliton production in 
particle collisions. 

In Fig. 22 we plot the semiclassical exponent at g'^N = 2: dashed lines represent the 
numerical data at different values of solid line is the result of extrapolation to 5p = 0. 
The results for Fn {E) at = 0 are also shown in Fig. 4b at different values of g'^N. One 
observes that the exponent decreases with energy approaching constant a.t E ‘^Ms- In [14] 
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we argued that this is the expected behavior for the collision-induced tunneling aX E ^ +oo. 

Note that numerical extrapolation to = 0 is harder to perform near the threshold 
Ecb = 2Ms because the asymptotic expansion of F/v in 6p has different form at smaller 
energies. However, the thin-wall arguments of Appendix E suggest near-threshold behavior 

En{E) ^ ci{N) + C 2 {N)\E - 2Ms |at E^ 2Ms , (4.14) 

where Ci and C 2 are unknown functions of N. We hnd that the numerical results in Figs. 22 
and 4b are consistent with the asymptotics (4.14) (dotted lines in both hgures). 


4.6 Two—particle processes 


Now, consider creation of the soliton-antisoliton pairs in the two-particle collisions {N = 2). 
These processes are natural from the viewpoint of collider physics but cannot be described by 
direct semiclassical methods. We compute their leading exponent E 2 {E) using the Rubakov- 
Son-Tinyakov conjecture (2.5), i.e. evaluating the limit g'^N —)■ 0 of the multiparticle expo¬ 
nent F/v(F). Recall that the semiclassical solutions develop a singularity in this limit. Thus, 
we cannot address them directly. In what follows we extrapolate the numerical results to 
g'^N —)■ 0 using an educated guess on the behavior of the multiparticle exponent. 



Figure 23: (a) Extrapolating the suppression exponent to g‘^N = 0 at g‘^E = 10 and 

6p = 0.06. (b) Two-particle exponent E 2 {E) at different dp (dashed lines) and the result of 
its linear extrapolation to = 0 (solid line). 


One hnds that vanishingly small initial particle number is achieved at 0 —>■ -(-cx): in 
this case the initial condition (2.2d) reduces to the Feynman one, and the initial state of 
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the process becomes semiclassically indistinguishable from the vacuum. Besides, since the 
combination e“® enters the semiclassical equations, one expects regular expansion 

g'^N = rii ■ e“® + 77-2 ■ e“^® + ... ^ 9 = - log{g‘^N) + 9o + 9i {g^N) + ... , (4.15) 


where n* and 6i are the energy-dependent Taylor coefficients. Now, we integrate the Legendre 
transform (2.6) at hxed energy. 


N 



0 


+00 



e{N) 


(4.16) 


where —>■ F 2 and g‘^9N —)■ 0 at g‘^N —)■ 0; besides, we integrated by parts in the second 
equality. Substituting the expected behavior (4.15) into Eq. (4.16), we obtain the asymptotic 
form 

F^iE) + g^NO ^ F^iE) - g^N + ■ {g^Nf + ... , (4.17) 

of the exponent at small g‘^N. This behavior, if conhrmed, automatically implies that the 
limit g‘^N —0 of the multiparticle exponent exists. 

In Fig. 23a we compare the numerical data for the quantity F^ + g‘^N6 (points) with the 
expectation (4.17) (dashed line). One observes a consistent £t involving two parameters, F 2 
and 9i. Changing the number of data points in the £t, we learn that the result for F 2 {E) is 
stable with relative precision of order 1%. 

The final numerical graphs of the two-particle exponent F 2 {E) are shown in Fig. 23b: 
dashed and solid lines are obtained at > 0 and dp —)■ 0, respectively. The last graph 
corresponding to soliton-antisoliton production is repeated in Fig. 4b. 


5 Concluding remarks 

In this paper we developed a new semiclassical method to calculate the probability of the 
topological soliton-antisoliton pair production in particle collisions. Our main idea was 
previously reported in [22] , here we presented the details, confirmations and generalizations 
of the method. 

In spite of all technicalities, the essence of our approach is simple: it is based on expecta¬ 
tion that the semiclassical solutions describing classically forbidden transitions continuously 
depend on parameters of the transitions.^^ With this idea in mind, we introduced a small 

general the continuity may be violated due to Stokes phenomenon [53, 43, 54]. We demonstrated, 
however, that the latter is not relevant for the processes considered in this paper. 
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external field £ coupled to the topological charges of the soliton and antisoliton. We started 
from the well-known solutions describing spontaneous Schwinger creation of the soliton pairs 
in the held £. Then we added colliding particles to the initial state of the process and ob¬ 
tained the solutions for the collision-induced Schwinger processes. Finally, we gradually 
increased the collision energy E above the threshold 2Ms of soliton pair production and 
switched off the held, T —)■ 0. In this way we continuously arrived at the solutions describing 
soliton-antisoliton production in particle collisions. 

We illustrated the above semiclassical method with the explicit numerical calculations 
in the (1 -|- l)-dimensional model of a scalar held. In particular, we demonstrated that 
the above-mentioned semiclassical solutions continuously depend on the parameters E and 
hp oc \£\ of the induced Schwinger process. We showed that these solutions reproduce the 
thin-wall results at low energies and correctly describe classical over-barrier transitions at 
E > 2Ms and high initial multiplicities. Finally, we computed the semiclassical suppres¬ 
sion exponents of the soliton-antisoliton pair production in the iV-particle and two-particle 
collisions, see Fig. 4b. 

We believe the semiclassical approach of this paper will be useful in other models / setups 
of particle and condensed matter physics. 

Acknowledgements. We thank V. Rubakov [55] for criticism and F. Bezrukov, D. Gor¬ 
bunov, M. Libanov, E. Nugaev, S. Troitsky, S. Sibiryakov for discussions. This work is 
supported by the RSCF grant 14-22-00161. Numerical calculations have been performed on 
the Computational cluster of the Theoretical division of INR RAS. We thank G. Rubtsov 
for supporting its stability. 

A Deriving the equations 

In this Appendix we review the semiclassical method of [28], see also [5, 37]. To this end we 
introduce path integral for the probability (2.1) and evaluate it in the saddle-point approx¬ 
imation. We obtain equations for the saddle-point conhguration x) and expression for 
the suppression exponent Fn{E). 

Semi-inclusive initial states in Eq. (2.1) have fixed energy E and multiplicity N. The 
projector Pe,n on the subspace of these states simplifies [28] in the coherent basis [56], 

p—ioo ( /‘I 

{a\pE,N\h) = - \^ET + N 9 + J dk'fkalbkj , . (A.l) 
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where |a), \b) are the eigenstates of the false-vacuum annihilation operators with eigen¬ 
values ttk, bk and normalization {a\b) = exp (J dkalbk). Integrals over T and 9 run along 
the imaginary axes. Expression (A.l) can be proven [37] by acting on the Fock states 
®Ii • • • coherent-state representation. We transform Eq. (A.l) into the configu¬ 

ration representation using 


{(j)i\a) = exp <{ - 


dk (-afcO.fc - ujk(j)i{k)(j)i{-k)/g‘^ 2y/^ak(j)i{k)/g) 


where (j)i{k) stands for the spatial Fourier transform of (j)i{x). Here and below the prefactors 
are omitted. We obtain the matrix elements 

p—iOQ 

(0i|^E,iv|0') = / . (a.2) 

J ioo 

denoting with 

^ = / 2g^{'f^- 1) ^ {(j)i{k)(j)i{-k) + 0'(fc)0'(-/c)) - 47fc0i(fc)0'(-/c)] , (A.3) 

the quadratic functional of (pi and 0'. 

Next, we write the probability (2.1) in the path integral form, 

Vn = J VPi VP[ VPf (0/|H(f/, ti)\PP) {cPi\PE,N\<t>[) U)\Pf) 

—ioo 

dTde 

XD 

In the hrst line of Eq. (A.4) we integrate over all initial and hnal-state conhgurations 0j, 0' 
and 0/ projecting onto the relevant subspace of initial states with Pe,n- In the second line 
the path integrals for the evolution operators U and f/I were introduced. We assume that 
the conhgurations 0 and 0' describe false vacuum decay: they are close to the false vacuum 
at f = fj —)■ —oo and contain expanding bubbles as t = t/ —)■ -|-oo. 

The integral (A.4) can be evaluated in the saddle-point approximation at large S and 
B i.e. at g^ -C 1, see Eqs. (1.3) and (A.3). In this case the dominant contribution to the 
integral comes from the vicinity of the saddle-point conhguration {ps(t, x), p'g(t, x), T, 9} 
which extremizes the integrand in Eq. (A.4). This means that the saddle-point helds 0^ and 
0(. satisfy the classical held equations 6S/6p = 6S/6p' = 0. Boundary conditions for the 
latter equations are obtained by varying the integrand with respect to 0* = 0(fj, x), 0' and 


VpiVp[Vpf [ Vp\^/ QiS[,j>]-iS[<f>']+2ET+Ne+Bl^i,4>'^ _ I'A.d) 
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i(pi{k) + Uk(pi{k) = 7fc yi^'iik) + Uk(p'i{k)j , 

(A.5) 

i(p'i{k) - Uk(p[{k) = 7fc (^i(pi{k) - Uk(pi{k)^ , 


0/ = 0/ , 0/ = 0/ , 

(A.6) 


where the dots denote the time derivatives coming from the variations of the action, e.g. 
6S/6(j)f = (pf/g^. Finally, differentiation with respect to T and 9 gives eqnations 


2E 


dB 

&T ’ 



(A.7) 


where B is dehned in Eq. (A.3). 

Two observations greatly simplify the semiclassical eqnations. First, cp' originates from 
the complex conjngate amplitnde snggesting 


(p'^{t,x) = [(ps{t, x)]* . (A.8) 

This Ansatz is compatible with Eqs. (A.5) if the saddle-point valnes of T and 9 are real, 
and we assnme that as well. Second, in the inhnite past the solntion (ps{t, x) describes free 
waves (2.2c) in the false vacunm. Hence, 

•#>.,.(*:) = [M + 9|-,| . (A.9) 

^/2uk ' ' 

At this point of calculation the initial time ti is real; it will be continued to complex values 
in Sec. 2. Substituting Eqs. (A.8), (A.9), into Eqs. (A.6), (A.5), (A.7), one obtains the 
boundary conditions (2.2b), (2. 2d) and expressions (2.3) relating (T, 9) to {E, N). 

We hnally substitute the saddle-point conhguration into the integrand of Eq. (A.4) and 
obtain the leading suppression exponent 

En{E) = -g\2ET + N9) + 2gHmS[(Ps] - g^B[(Ps,^, (P*J . (A.IO) 

Rewriting the last term with help of Eq. (A.9) and the boundary condition (2. 2d), we arrive 
at the standard expression (2.4). 

We hnish derivation with two remarks. First, the solution (ps and all equations will be 
considered along the complex-time contour in Fig. 5 corresponding to tunneling. In this 
case the initial time ti is taken complex because the exponent (2.4) does not depend on it 
anyway. Second, the functional (A.IO) depends on E, N explicitly and via the saddle-point 
conhguration {(pg, T, 9}. However, since we have already extremized this functional with 
respect to (p, T and 9, its E and A^-derivatives are given by Eqs. (2.6). 


40 




B Lattice formulation 


In this Appendix we list the discretized action and finite-difference seniiclassical equa¬ 
tions (2.2), see [11, 12] for similar approaches. 

Our lattice {tj, Xi] is defined in Sec. 3.2. Performing the substitutions (3.2) and (3.3) in 
the action (1.3), we obtain. 


9^S = 


Nt 

N^-1 

E 

i=0 

j=-i 


Ax, 




At, 


Nt+l 

N^-2 

E 

i=0 

j=0 








Ax 


Nt+l 

N^-l 

2 Atj Axi V (0j^, 

i=0 

j=0 


(B.l) 


Recall that the lattice field equations are the derivatives of this expression at the internal 
points 0 < j < iVt of the time contour. Assuming for simplicity 1 < i < — 2, we obtain. 


jr ^ 0i+i 
-rj,i — 








h', *+l 


+ 


h'q-1 




Atj Atj 


Atj_iAtj 




+ R'(0,-,) = O, (B.2) 


where V is the derivative of the scalar potential. Equations at the spatial boundaries, i.e. 
at i = 0 and i = — 1, are derived in similar manner. 

We assume that the evolution is linear in the beginning of the process. In this case the 
potential reduces to R ~ m?_(0 — 0_)^/2. Substituting it into the above action, we obtain 
Eq. (3.6) with discrete time and 





^i, i' 


Si+l,i' + 

Ax^J AxiAx,! 


We directly compute the real eigenvectors and eigenvalues of this symmetric matrix 
by QR decomposition. The coefficients /„ and gn of the linear solution (3.7) are then given 
by Eqs. (3.8), where the time derivative (3.9) simplifies. 


5 ^ dtA 



'ipijto) - 'ipijt-i) 
At_i 


At_iujl'ijji{t.i) 


(B.3) 


Using explicit expressions for fn and Qn, we rewrite the initial condition fn = ^ in the 
matrix form. 


X^d^^Kl-e ^)uJnRe'ilJi- {1 + e = 0 , 

i 

+ + (1 - e"^)Re9i'0i}j_^ = 0 , 

i 
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where the discrete operator dt is dehned in Eq. (B.3). One reads off the matrices Mr and 
Mj in Eq. (3.10) from the above equations. 

We directly discretize the nonlinear energy (3.11), 


N^-l 


g‘^Ej+i/2 = ^ Ax, 




2=1 


'N^-2 

i=0 


(0jq+i Y j 

2 Ax 




- -h ^ Axi V{(j)j^i) + (j —t j + 1) z • (B.4) 


i=0 


It conserves, i.e. does not depend on j up to O(At^) numerical errors. Somewhat different 
discretization procedure is natural for the energy of the linear system (3.6) at t = t_i. We 
start from the expression 


N^-l 


= 




i=0 




N^-1 

E 

2 , 2^=0 


V^ 2 -^ 2 , 


(B.5) 


in continuous time and discrete space. Substituting the solution (3.7), we obtain the hrst of 
Eqs. (3.12). The respective formula for the initial multiplicity is then easily deduced from 
Eqs. (2.3). 

We hnish this Appendix with the lattice expression for the suppression exponent^'^ 
Fr,{E) = -g\2ET + NO) + 2^721111 + Im Ax, (0_i,, + - 20_) , (B. 6 ) 

i=0 

where the action is given by Eq. (B.l). 


C Tests of the numerical procedure 

We subdued the lattice solutions to a number of consistency tests which support the numer¬ 
ical methods of Sec. 3 and allow us to estimate the numerical errors. 

For a start, we checked sensitivity of the results to the spatial cutoff L^. To this end we 
increased the cutoff from = 7 to 14 at a hxed lattice spacing Ax = L^/{Nx — 1 ). The 
integral quantities E, N and F/v stayed independent of up to relative errors of order 10“^. 

last integral in this expression is 0{At) and 0{Ax‘^) accurate. In practice the first-order correction 
in At is negligible because |At| ^ Ax. All other lattice expressions in our study are second-order both in 
space and time. 
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Figure 24: (a) Full and linear energies of the semiclassical solution with {E, N) ^ 
(9.0, 2.8)/g‘^ at Ret < 0. (b) The exponent F/v as a function of N at different lattice 
spacings. The space-time box in both figures is = 7, Re (tNt+i — t-i) ~ 15. 


Another source of numerical errors is related to the hnite extent t_i < t < tAr^+i of the 
temporal lattice. Recall that the semiclassical solutions should describe free waves in the false 
vacuum at t ~ t_i. This property is conceptually important, it was used in the derivation of 
the initial condition (2.2d). We estimate the effect of nonlinear interactions in the beginning 
of the process by evaluating the energy of the linearized system given by Eqs. (3.12), (3.8) 
at different time sites tj and comparing it with the full conserved energy (3.11), see Fig. 24a. 
The graphs stay close at t ~ —8.5, separating in the nonlinear region t > —3. The 

relative error due to nonlinear interactions in the beginning of the process is estimated as 
|Ffuii — Fiinearl/Ffuii < 10“^. One concludes that the semiclassical evolution at t ~ is, 
indeed, free thanks to the clever choice of the scalar potential in Sec. 3.1. 

Let us turn to discretization effects which should be O(Af^) and 0{Ax‘^) small in the 
second-order finite-difference scheme. First of all, one observes that conservation of the full 
energy in Fig. 24a is violated at the level of AFfun/Ffuu < 10“^. This is the effect of time 
discretization because energy conservation is restored at At —)■ 0 and finite Ax. Second, we 
directly estimate the finite-difference errors comparing numerical results at different lattices, 
see Figs. 24b and 25. Recall that our reference-point lattices are different at energies below 
and above Ecb'- Nt x = 3200 x 150 and 11000 x 500, respectively. In the former case the 
relative errors in F, N and Fjv stay below 1% reaching maximum at the smallest N. Fine 
lattice resolution at high energies gives errors well below 1% at F ~ E^t ~ 6.1 and high N. 
The errors grow, however, to 1% at F > Id/gf^ and/or small multiplicities. In what follows 
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Figure 25: The integral quantities E, N and F/v as functions on the (a) time and (b) spatial 
lattice spacings at (T, 9) = (0.026, 0.5) and (F, N) k. (7.8, 2.b)/g‘^] E > Ecb- The dashed 
lines are the linear hts of the data points. We use the space-time box from Fig. 24. 
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Figure 26: The discrepancy in the initial condition arg /2 = arg 5 f 2 at T = 0.026 and differ¬ 
ent 6] E > Ecb- The standard lattice size and the space-time box from Fig. 24 are used. 

we exclude the results with g'^E > 14 and g‘^N < 1 due to improperly high discretization 
errors. 
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Figure 27: Partial derivatives of the suppression exponent with respect to the (a) energy E 
and (b) initial particle number N. The graphs are plotted along the lines (a) T = 0.026 and 
(b) e = 0.5. 


To conclude, we keep the hnite-difference effects below the relative level of 1%. Figure 25 
shows that most of the integral quantities linearly depend on oc N~‘^ and oc 
like they should in the second-order discretization scheme. The only exception is the energy 
E (lower panel in Fig. 25b) which receives small nonpolynomial correction from adiabatic 
high-frequency waves in the solution. This effect is negligible in the region E < which 

we consider (see, however, the study of the high-energy solutions in [14]). 

In Sec. 3.3 we traded the initial condition arg /„q = arg^f^g at a given hq for the artihcial 
constraint hxing the time-translation invariance. We argued that the omitted condition will 
be automatically satished once the other lattice equations are solved correctly. In Fig. 26 we 
demonstrate that this is, indeed, the case: at no = 2 the related numerical error is smaller 
than 10“^. 

A good piece of our qualitative and quantitative results is based on the Legendre trans¬ 
forms (2.6) which should hold with acceptable numerical precision. In Fig. 27 we plot the 
partial derivatives {OeEj^, dNEN)/g‘^ of the suppression exponent and compare them to —2T 
and —9 (dashed horizontal lines). One observes that the Legendre transforms hold with ab¬ 
solute precision AT, A9 < 10“^. 

To summarize, the relative hnite-volume and discretization effects in our solutions are 
smaller than 10“^ and 10“^, respectively. The other numerical errors are vanishingly small. 
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D Solving the linear system 


Each iteration of the Newton-Raphson method involves solution to the system of {Nt + 3) x 

> 10® complex linear equations (3.16). This is the major CPU time-consuming part of 
the numerical procedure. We decrease its computational cost using the sparse structure of 
the lattice held equations. 

In what follows we suppress the spatial indices working with the V 2 .-dimensional vector 
held (j)j = ■ ■ ■, Recall that the linear system (3.16) is obtained by substituting 

(pj = 0^°^ + 6(pj into the lattice equations and expanding them to the linear order in Scpj. In 
particular, the discrete held equations (B.2) take the form, 

Cj • S(j)j = Lj • T Rj • — Rj , (D.l) 

where Rj is their left-hand side at 0 = 0^®\ while Cj, Lj and Rj are the x coefficient 
matrices which can be deduced from Eqs. (B.2), e.g. (Cj)j^j/ = dJ-'j^i/d(l)j^i'\^(^gy Note that 
Eq. (D.l) is sparse: it relates (j0j-i, Scpj and Scpjyi. Besides, the matrices Lj and Rj are 
diagonal, while Cj is three-diagonal. We will use both these facts while solving the linear 
system. 

Our “stable” algorithm [31] eliminates equations from the set (D.l). Namely, suppose 
we express Scpj from the j-th equation, 

6(j)j = C-^ {Lj ■ (50,_i + R, ■ 50,+1 - Rj) , (D.2) 

and substitute it into the neighbouring {j =p l)-th equations. The latter will keep the form 
(D.l) albeit with new coefficient matrices. In particular, the {j — l)-th equation will relate 
50,■_ 2 , 50,-1 and 50,+i with 

= Q_i - R,_i C.^ L, , = L,_i , 

R'j-i = Rj-i C~^ Rj , R'j_^ = Rj-i + Rj-i Cj^ Rj . (D.3) 

Repeatedly using Eqs. (D.3), we eliminate all held equations except for the very hrst and 
last ones at j = 0, W- To make this procedure more stable, we hrst apply (D.3) to the 
equations with odd j, then eliminate odd equations from the remaining set, etc. Once we 
are left with the equations at j = 0 and W relating 50_i, 50o, 50Arj, 507 Vt+i, we add the 
linearized boundary conditions and solve the resulting linear system by the direct method 
of LU decomposition. The complete solution {50,} is restored from the boundary values of 
50 and Eqs. (D.2). 
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The above elimination process is apparently ineffective. Indeed, the substitutions (D.3) 
do not preserve the sparse form of the coefficient matrices thus requiring general matrix 
multiplications and inversions at the second stage of the elimination process. One concludes 
that cNtN^ operations with c ~ 1 are needed for obtaining the solution. Note, however, 
that the disadvantages of this “slow” algorithm are compensated by its exceptional stability 
properties. We exploit it at < Ecb where the faster procedure accumulates round-off errors 
and causes divergence of the Newton-Raphson iterations. 

Our “fast” method benefits from the sparse form of the coefficient matrices in Eq. (D.l). 
It is based on Cauchy problem for the lattice field equation. Consider the x (2W + 1) 
matrix satisfying the analog of Eq. (D.l), 

~ ~ ^i) > 

where the last matrix term in the right-hand side contains the vector Ej in the last column; 
the bold 0 and 1 denote zero and unit W x W matrices. We solve the Cauchy problem for 
the “propagator” H with the initial conditions 

S _1 = (1, 0, 0) , So = (0, 1, 0) . (D.5) 


This procedure involves multiplications of Hj by the diagonal and three-diagonal matrices 
Lj, and Cj, i.e. cNtN^ operations, c ~ 1. 

Given the propagator S, we introduce a convenient representation of the solution. 


6(j)j 


A0-l\ 
S(j)o , 

\ 1 / 


(D.6) 


where the {2Nx + l)-vector in the right-hand side is composed from the boundary fields 
6(p-i and S(j)Q. One can check that the vector (D.6) satisfies Eq. (D.l). Equation (D.6) 
relates, in particular, to S(f)-i and 6(f)o. Adding these two relations to the set of 

boundary conditions, we obtain a closed linear system for 6(j)o, S(j)Nt and 6(j)Nt+i- We 

solve the latter using the LU decomposition method. Given and we Gauchy-evolve 
Eq. (D.l) determining the solution {S(j)j}. 

In the “fast” method we arrive at the solution in cN^N^ operations, a factor of faster 
than in the “slow” algorithm. The AC and AC scalings of the GPU times in our “slow” and 
“fast” codes are illustrated in Figs. 28a,b. 

Let us point at the reason behind the poor stability properties of the “fast” algorithm at 
low energies. We argued in Sec. 4 that the distinctive property of the semiclassical solutions 
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Figure 28: Scaling with the lattice size of the CPU time required for: (a) elimination of 
the equations in the “slow” algorithm, (b) computation of the propagator S in the “fast” 
one. Eight processors are used. Regard drastically different lattice ranges in Figs, (a), (b): 
Nt = 3200, = 100 ^ 250 and Nt = 11000, = 250 ^ 2000, respectively. 


at E < Ecb is long periods T of their Euclidean evolutions. Linear perturbations 6(j)j grow 
exponentially in Euclidean spacetime magnifying the initial round-off errors. They cannot be 
correctly evolved within the “fast” Cauchy approach as opposed to the homogeneous “slow” 
algorithm. As a consequence, we exploit the “fast” procedure only for almost-Minkowskian 
solutions at E > Ecb- 

Both our algorithms can be performed in parallel. To this end one divides equations (D.l) 
into the subsets with index ranges j = 0 ... N^, ... 2Nk, etc., and performs computations 

in every subset at a separate processor. At the second stage of the algorithm one combines 
the data from all processors. For example, in the parallel version of the “fast” algorithm the 
propagators in all subsets are computed, and the original propagator S is restored as 
their product. In Fig. 29 we demonstrate that the elapsed real time for obtaining S scales 
as c/Nproc with the number of processors. 


E Thin—wall approximation 

The semiclassical solutions describing false vacuum decay can be found analytically if the 
sizes of their true vacuum bubbles are parametrically large compared to the widths of the 
bubble walls. We will see that this thin-wall regime [3, 33, 10, 36] occurs at —)■ 0 
and E < Ecb- 
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Figure 29: The time for computing the propagator S at Nproc = 2d-ll processors; Nt x = 
11000 X 250. The data points scale as c/Nproc (dashed line). 


Consider some Euclidean held conhguration with a large true vacuum bubble, see Fig. 30a. 
At a crude level it can be characterized with the positions of the bubble walls xi{t) and 
3^2 ('r), where r = if is a Euclidean time. The corresponding Euclidean action receives large 
contributions from the constant energy density {—Sp/g"^) inside the bubble and tension of 
its walls, cf. Eq. (1.3), 


Se ~ MsL - 5pA/g‘^ 



M* 


'sj\ x\ + Mg\ x 


6p (X2 - Xi)/g‘^ 


(E.l) 


Here A is the spacetime area occupied by the true vacuum, L is the length of its boundaries; 
in the second equality we expressed A and L in terms xi{t), X 2 {t) denoting the r-derivatives 
with dots. Anticipating the limit 6p —)■ 0, we identihed the wall tension with the soliton mass 
Ms- In what follows we treat the action (E.l) as a functional of Xi(r) and X 2 {t). 

Now, let us construct the periodic instantons, i.e. solutions to the Euclidean held equa¬ 
tions with the boundary conditions (4.1). One notes that these are the extrema of the 
Euclidean action in the interval r G [—T, 0]. Indeed, the Neumann conditions (4.1) mean 
that the action is stationary with respect to the boundary values of the helds at r = —T 
and 0. Working in the thin-wall approximation, we extremize the action (E.l) by varying 
Xi(r) and X 2 {t). One can draw some intuition here by noting that the latter action is similar 
to the static energy of a soap bubble in two space dimensions. Its extremum is achieved if 
the bubble walls Xi and X 2 form circular arcs with hxed curvature radius 


Rb = 


9‘^Ms 

6p 
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(b) 


(c) 


Figure 30: Thin-wall configurations in Euclidean spacetime (r, x): general configuration (a) 
and the stationary solutions with (b) T < Rb and (c) T > Rb. Pink/gray and white regions 
are hlled with the true and false vacua, respectively. 

Extremization with respect to the boundary values of xi and X 2 gives conditions xi = X 2 = 0 
at r = —T and 0. The only exception appears if the arcs Xi and X 2 meet at one point as in 
Fig. 30b; then the boundary condition changes to = —X 2 - 

The stationary thin-wall conhgurations ai T < Rb and T > Rb are shown in Figs. 30b 
and 30c, respectively. Substituting them into the action (E.l), one obtains. 


Sb 


g^Mg \ 2Q; + sin2Q; at T = sin a < , 

25p I TT at T > Rb , 


(E.2) 


where a is the half-angle between the arcs in Fig. 30b. One obtains Eq. (4.6) for the 
Minkowski action S = iSE in terms of the half-period T. 

Let us hnd out when the thin-wall approximation is trustworthy. One naively expects 
this regime to be solid at small 6p because all bubble sizes are proportional to Rb oc l/6p. 
However, the distance between the bubble walls in Fig. 30b vanishes at a —?• 0 as a^Rb. 
Comparing it with the typical wall width mT^, one obtains the correct condition 


a > (m+Rb) , 


(E.3) 


for the thin-wall approximation. This inequality breaks down when the solution approaches 
the critical bubble. Indeed, consider the thin-wall energy 


E = 2Ms cos a , 


(E.4) 
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which is easily deduced from Eq. (E.l) and Fig. 30b. It tends to 2Ms as a, T —)■ 0. Thus, at 
small Sp the thin-wall approximation is trustworthy below the barrier height E^b but breaks 
at E > Ecb- 

Next, we consider transitions from the semi-inclusive initial states with fixed energy E 
and multiplicity N. Note that the periodic instantons describe particular processes of this 
kind with N = Npi{E), see Fig. 11. Their in-states are obtained by continuing the solutions 
from T = it = —T to the initial part AB of the time contour in Fig. 5a. One notices [36], 
however, that the r = —T sections of the thin-wall configurations in Figs. 30b, c do not 
depend on 6p. This suggests that the initial states of the periodic instantons and their 
matrix elements with the states of smaller multiplicities are also independent of Sp. Thus, 
decreasing the initial multiplicity from N = Npj to N = 2 one at best gets 0{Sp^) correction 
to the exponent Epfi^E). We conclude^® that the leading 1/Sp part of En{E) does not depend 
on N. Expressing the action (E.2) in terms of energy (E.4) and substituting it into Eq. (2.4), 
one obtains the thin-wall result (4.7), (4.8) for the suppression exponent at F < Ecb- 

Finally, let us guess the behavior of En{E) near the point E ^ 2Ms where the asymptotic 
expansion in Sp breaks down. Inequality (E.3) shows that the thin-wall approximation is 
valid at 

where we used Eq. (E.4) and m+ ~ g^Ms- One expects that the terms of the thin-wall 
expansion (4.7) become comparable at the boundary of this interval. In particular. 


E. 


N,0 


Sp 3Sp 


1/2 


\2Ms - E 


2,12 


where the explicit form (4.8) was used. We express Sp from Eq. (E.5), and obtain Ejsifl ~ 
const ■ \E — 2Ms\^^‘^. Assuming some regular contribution besides this singular term, we 
arrive at Eq. (4.14). 
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